Data Analysis: Acute Changes in CEWL as a Function of Temperature

Calvin Davis, Savannah Weaver

2022

Packages

if (!require("tidyverse")) install.packages("tidyverse") 
library("tidyverse") #tidyr, ggplot, dplyr, etc
if (!require("ggpubr")) install.packages("ggpubr") 
library("ggpubr")
if (!require("lme4")) install.packages("lme4") 
library("lme4") # mixed-effect models
if (!require("lmerTest")) install.packages("lmerTest") 
library("lmerTest") # p-values
if (!require("car")) install.packages("car") 
library("car") # test for VIF
if (!require("broom.mixed")) install.packages("broom.mixed")
library("broom.mixed") # model export
if (!require("emmeans")) install.packages('emmeans')
library('emmeans')
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # model export
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown format

Background

This code walks through the statistical models and figures for a study on the rate of change of cutaneous evaporative water loss (CEWL) in response to temperature change in Western Fence Lizards (Sceloporus occidentalis) (published in the Journal __ in 2023). Data was collected and analyzed by Calvin Davis and Savannah Weaver, under the advising of Dr. Emily Taylor at California Polytechnic State University.

Load Data

Get the RDS files created in the data wrangling Rmd:

temp_time_series <- read_rds("Data/temp_time_series.RDS")
CEWL_time_series <- read_rds("Data/CEWL_time_series.RDS")

Stats & Models

Temp Correlation

Check the correlation between dorsal and cloacal temps:

temp_corr <- lm(data = temp_time_series, 
                calibrated_dorsal_temp ~ calibrated_cloacal_temp)
summary(temp_corr)
## 
## Call:
## lm(formula = calibrated_dorsal_temp ~ calibrated_cloacal_temp, 
##     data = temp_time_series)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0389 -0.7784  0.0075  0.5911  5.3351 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             6.966310   0.049820   139.8   <2e-16 ***
## calibrated_cloacal_temp 0.747869   0.001899   393.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.388 on 25249 degrees of freedom
## Multiple R-squared:   0.86,  Adjusted R-squared:  0.8599 
## F-statistic: 1.55e+05 on 1 and 25249 DF,  p-value: < 2.2e-16
#write.csv(broom.mixed::tidy(temp_corr), 
 #         "./Results_Statistics/body_temp_correlation.csv")

They are actually pretty different, so run parallel models of CEWL, one with dorsal temp, and one with cloacal temp.

CEWL (raw) ~ temp LMM

First, test linear effects of temperature:

CEWL_LMM_02a <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             calibrated_cloacal_temp *
                             treatment +
                             (1|date/lizard_ID))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
summary(CEWL_LMM_02a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## CEWL_g_m2h ~ calibrated_cloacal_temp * treatment + (1 | date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 87927.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1047 -0.5895 -0.0371  0.4647  8.7021 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 10.066   3.173   
##  date           (Intercept) 22.540   4.748   
##  Residual                    1.882   1.372   
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                           Estimate Std. Error t value
## (Intercept)                               51.20985    2.56934   19.93
## calibrated_cloacal_temp                   -1.31478    0.03945  -33.33
## treatmentCooling                         -51.86569    1.71910  -30.17
## treatmentHeating                         -39.66638    1.78910  -22.17
## calibrated_cloacal_temp:treatmentCooling   1.62627    0.03957   41.09
## calibrated_cloacal_temp:treatmentHeating   1.54114    0.03957   38.95
## 
## Correlation of Fixed Effects:
##             (Intr) clbr__ trtmnC trtmnH cl__:C
## clbrtd_clc_ -0.421                            
## tretmntClng -0.472  0.628                     
## tretmntHtng -0.459  0.605  0.669              
## clbrtd_c_:C  0.419 -0.997 -0.629 -0.603       
## clbrtd_c_:H  0.419 -0.997 -0.626 -0.607  0.994
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
vif(CEWL_LMM_02a)
##                                         GVIF Df GVIF^(1/(2*Df))
## calibrated_cloacal_temp           309.047830  1       17.579756
## treatment                           1.853502  2        1.166805
## calibrated_cloacal_temp:treatment 311.665946  2        4.201674
CEWL_LMM_02b <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             calibrated_dorsal_temp *
                             treatment +
                             (1|date/lizard_ID))
summary(CEWL_LMM_02b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ calibrated_dorsal_temp * treatment + (1 | date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 88772.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9809 -0.5887 -0.0570  0.4961  8.9013 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 10.316   3.212   
##  date           (Intercept) 21.015   4.584   
##  Residual                    1.946   1.395   
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                              51.27388    2.54285   20.16
## calibrated_dorsal_temp                   -1.30362    0.04142  -31.48
## treatmentCooling                        -53.18221    1.77431  -29.97
## treatmentHeating                        -40.88244    1.84485  -22.16
## calibrated_dorsal_temp:treatmentCooling   1.64631    0.04157   39.60
## calibrated_dorsal_temp:treatmentHeating   1.57709    0.04164   37.88
## 
## Correlation of Fixed Effects:
##             (Intr) clbr__ trtmnC trtmnH cl__:C
## clbrtd_drs_ -0.451                            
## tretmntClng -0.499  0.645                     
## tretmntHtng -0.486  0.623  0.680              
## clbrtd_d_:C  0.449 -0.996 -0.647 -0.620       
## clbrtd_d_:H  0.449 -0.995 -0.642 -0.625  0.991
vif(CEWL_LMM_02b)
##                                        GVIF Df GVIF^(1/(2*Df))
## calibrated_dorsal_temp           231.169627  1       15.204263
## treatment                          1.944891  2        1.180929
## calibrated_dorsal_temp:treatment 234.207111  2        3.912011

These are REALLY good models, but based on the figures below, the relationship isn’t actually linear… so test the inclusion of polynomial factors:

CEWL_LMM_03a <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             poly(calibrated_cloacal_temp, 2)*treatment +
                             (1|date/lizard_ID))
summary(CEWL_LMM_03a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 81611.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6892 -0.6401 -0.0881  0.5378  7.8576 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 11.773   3.431   
##  date           (Intercept) 23.794   4.878   
##  Residual                    1.468   1.212   
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                                      Estimate Std. Error
## (Intercept)                                            0.6557     2.4503
## poly(calibrated_cloacal_temp, 2)1                   1422.2807    62.9113
## poly(calibrated_cloacal_temp, 2)2                  -3419.1801    82.4963
## treatmentCooling                                       6.6956     1.5010
## treatmentHeating                                      16.7664     1.5897
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling -1162.0238    62.9471
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling  3565.6587    82.5285
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating -1223.6577    62.9510
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating  3338.2055    82.5256
##                                                    t value
## (Intercept)                                          0.268
## poly(calibrated_cloacal_temp, 2)1                   22.608
## poly(calibrated_cloacal_temp, 2)2                  -41.446
## treatmentCooling                                     4.461
## treatmentHeating                                    10.547
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling -18.460
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling  43.205
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating -19.438
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating  40.451
## 
## Correlation of Fixed Effects:
##             (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.159                                                        
## ply(c__,2)2  0.164 -0.914                                                 
## tretmntClng -0.336  0.260    -0.268                                       
## tretmntHtng -0.324  0.244    -0.251     0.505                             
## pl(__,2)1:C  0.159 -0.999     0.914    -0.259 -0.244                      
## pl(__,2)2:C -0.164  0.914    -1.000     0.268  0.251 -0.913               
## pl(__,2)1:H  0.159 -0.999     0.914    -0.259 -0.244  0.999     -0.913    
## pl(__,2)2:H -0.164  0.914    -1.000     0.268  0.251 -0.914      0.999    
##             p(__,2)1:H
## ply(c__,2)1           
## ply(c__,2)2           
## tretmntClng           
## tretmntHtng           
## pl(__,2)1:C           
## pl(__,2)2:C           
## pl(__,2)1:H           
## pl(__,2)2:H -0.914
vif(CEWL_LMM_03a)
##                                                    GVIF Df GVIF^(1/(2*Df))
## poly(calibrated_cloacal_temp, 2)           9.312576e+05  2       31.064721
## treatment                                  1.100911e+00  2        1.024326
## poly(calibrated_cloacal_temp, 2):treatment 9.312225e+05  4        5.573547
anova(CEWL_LMM_03a, CEWL_LMM_02a)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_02a: CEWL_g_m2h ~ calibrated_cloacal_temp * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_02a    9 87929 88002 -43956    87911                         
## CEWL_LMM_03a   12 81676 81774 -40826    81652 6259.3  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
CEWL_LMM_03b <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             poly(calibrated_dorsal_temp, 2)*treatment +
                             (1|date/lizard_ID))
summary(CEWL_LMM_03b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 83028.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7628 -0.6345 -0.0833  0.5160  8.3769 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept)  9.996   3.162   
##  date           (Intercept) 21.953   4.685   
##  Residual                    1.554   1.246   
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                                    Estimate Std. Error t value
## (Intercept)                                           3.248      2.336   1.390
## poly(calibrated_dorsal_temp, 2)1                   1534.609     67.430  22.759
## poly(calibrated_dorsal_temp, 2)2                  -2711.185     75.106 -36.098
## treatmentCooling                                      3.529      1.386   2.546
## treatmentHeating                                     14.341      1.468   9.766
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling -1350.765     67.457 -20.024
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling  2852.577     75.134  37.967
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating -1365.961     67.469 -20.246
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating  2666.790     75.141  35.490
## 
## Correlation of Fixed Effects:
##             (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.162                                                        
## ply(c__,2)2  0.163 -0.946                                                 
## tretmntClng -0.327  0.273    -0.276                                       
## tretmntHtng -0.315  0.257    -0.260     0.507                             
## pl(__,2)1:C  0.162 -1.000     0.946    -0.273 -0.257                      
## pl(__,2)2:C -0.163  0.946    -1.000     0.276  0.259 -0.946               
## pl(__,2)1:H  0.162 -0.999     0.946    -0.273 -0.257  0.999     -0.945    
## pl(__,2)2:H -0.163  0.946    -1.000     0.276  0.259 -0.945      0.999    
##             p(__,2)1:H
## ply(c__,2)1           
## ply(c__,2)2           
## tretmntClng           
## tretmntHtng           
## pl(__,2)1:C           
## pl(__,2)2:C           
## pl(__,2)1:H           
## pl(__,2)2:H -0.945
vif(CEWL_LMM_03b)
##                                                   GVIF Df GVIF^(1/(2*Df))
## poly(calibrated_dorsal_temp, 2)           5.528187e+05  2       27.267523
## treatment                                 1.107683e+00  2        1.025897
## poly(calibrated_dorsal_temp, 2):treatment 5.527949e+05  4        5.221803
anova(CEWL_LMM_03b, CEWL_LMM_02b)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_02b: CEWL_g_m2h ~ calibrated_dorsal_temp * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_02b    9 88775 88848 -44379    88757                         
## CEWL_LMM_03b   12 83092 83190 -41534    83068 5689.1  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The poly models have significantly better explanatory power than the linear models.

Also try log transformation to make sure a polynomial effect is best:

CEWL_LMM_04a <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             (log(calibrated_cloacal_temp))*treatment +
                             (1|date/lizard_ID))
summary(CEWL_LMM_04a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ (log(calibrated_cloacal_temp)) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 88908.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1896 -0.5896 -0.0624  0.4489  8.9787 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept)  9.877   3.143   
##  date           (Intercept) 22.480   4.741   
##  Residual                    1.958   1.399   
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                               Estimate Std. Error t value
## (Intercept)                                    129.257      4.312   29.97
## log(calibrated_cloacal_temp)                   -34.468      1.097  -31.41
## treatmentCooling                              -143.864      3.871  -37.16
## treatmentHeating                              -131.365      3.906  -33.63
## log(calibrated_cloacal_temp):treatmentCooling   41.258      1.100   37.50
## log(calibrated_cloacal_temp):treatmentHeating   40.503      1.101   36.80
## 
## Correlation of Fixed Effects:
##             (Intr) lg(__) trtmnC trtmnH l(__):C
## lg(clbrt__) -0.842                             
## tretmntClng -0.844  0.938                      
## tretmntHtng -0.838  0.930  0.930               
## lg(clb__):C  0.840 -0.998 -0.940 -0.928        
## lg(clb__):H  0.840 -0.997 -0.935 -0.933  0.995
vif(CEWL_LMM_04a)
##                                             GVIF Df GVIF^(1/(2*Df))
## log(calibrated_cloacal_temp)           382.13042  1       19.548156
## treatment                               11.39233  2        1.837186
## log(calibrated_cloacal_temp):treatment 416.69458  2        4.518086
anova(CEWL_LMM_04a, CEWL_LMM_02a)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_04a: CEWL_g_m2h ~ (log(calibrated_cloacal_temp)) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_02a: CEWL_g_m2h ~ calibrated_cloacal_temp * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)
## CEWL_LMM_04a    9 88930 89003 -44456    88912                     
## CEWL_LMM_02a    9 87929 88002 -43956    87911 1000.4  0
anova(CEWL_LMM_04a, CEWL_LMM_03a)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_04a: CEWL_g_m2h ~ (log(calibrated_cloacal_temp)) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_04a    9 88930 89003 -44456    88912                         
## CEWL_LMM_03a   12 81676 81774 -40826    81652 7259.7  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
CEWL_LMM_04b <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             log(calibrated_dorsal_temp)*treatment +
                             treatment +
                             (1|date/lizard_ID))
summary(CEWL_LMM_04b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ log(calibrated_dorsal_temp) * treatment + treatment +  
##     (1 | date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 89739.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0661 -0.5838 -0.0705  0.4871  9.0055 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 10.302   3.210   
##  date           (Intercept) 21.092   4.593   
##  Residual                    2.024   1.423   
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                              Estimate Std. Error t value
## (Intercept)                                   131.771      4.497   29.30
## log(calibrated_dorsal_temp)                   -35.114      1.169  -30.05
## treatmentCooling                             -151.142      4.117  -36.71
## treatmentHeating                             -137.198      4.155  -33.02
## log(calibrated_dorsal_temp):treatmentCooling   43.244      1.172   36.89
## log(calibrated_dorsal_temp):treatmentHeating   42.178      1.174   35.94
## 
## Correlation of Fixed Effects:
##             (Intr) lg(__) trtmnC trtmnH l(__):C
## lg(clbrt__) -0.863                             
## tretmntClng -0.863  0.942                      
## tretmntHtng -0.857  0.934  0.933               
## lg(clb__):C  0.860 -0.997 -0.945 -0.931        
## lg(clb__):H  0.859 -0.996 -0.938 -0.938  0.993
vif(CEWL_LMM_04b)
##                                            GVIF Df GVIF^(1/(2*Df))
## log(calibrated_dorsal_temp)           275.88918  1       16.609912
## treatment                              12.81722  2        1.892119
## log(calibrated_dorsal_temp):treatment 315.47270  2        4.214446
anova(CEWL_LMM_04b, CEWL_LMM_02b)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_04b: CEWL_g_m2h ~ log(calibrated_dorsal_temp) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_02b: CEWL_g_m2h ~ calibrated_dorsal_temp * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)
## CEWL_LMM_04b    9 89762 89835 -44872    89744                     
## CEWL_LMM_02b    9 88775 88848 -44379    88757 987.18  0
anova(CEWL_LMM_04b, CEWL_LMM_03b)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_04b: CEWL_g_m2h ~ log(calibrated_dorsal_temp) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_04b    9 89762 89835 -44872    89744                         
## CEWL_LMM_03b   12 83092 83190 -41534    83068 6676.3  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The logarithmic models are NOT better than the linear or polynomial ones.

VIFs are still close to 1 with the polynomial included, and all variables seem significant, so check assumptions:

plot(CEWL_LMM_03a)

plot(CEWL_LMM_03b)

They both have a really weird spike, but otherwise there’s no fanning or different pattern.

Re-run with p-values:

CEWL_LMM_besta <- lmerTest::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             poly(calibrated_cloacal_temp, 2) * treatment +
                             (1|date/lizard_ID))
summary(CEWL_LMM_besta)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 81611.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6892 -0.6401 -0.0881  0.5378  7.8576 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 11.773   3.431   
##  date           (Intercept) 23.794   4.878   
##  Residual                    1.468   1.212   
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                                      Estimate Std. Error
## (Intercept)                                            0.6557     2.4503
## poly(calibrated_cloacal_temp, 2)1                   1422.2807    62.9113
## poly(calibrated_cloacal_temp, 2)2                  -3419.1801    82.4963
## treatmentCooling                                       6.6956     1.5010
## treatmentHeating                                      16.7664     1.5897
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling -1162.0238    62.9471
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling  3565.6587    82.5285
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating -1223.6577    62.9510
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating  3338.2055    82.5256
##                                                            df t value Pr(>|t|)
## (Intercept)                                            5.4550   0.268 0.798840
## poly(calibrated_cloacal_temp, 2)1                  25235.9702  22.608  < 2e-16
## poly(calibrated_cloacal_temp, 2)2                  25238.1540 -41.446  < 2e-16
## treatmentCooling                                      30.3218   4.461 0.000104
## treatmentHeating                                      29.8965  10.547 1.36e-11
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling 25235.9622 -18.460  < 2e-16
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling 25238.1521  43.205  < 2e-16
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating 25235.9530 -19.438  < 2e-16
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating 25238.1526  40.451  < 2e-16
##                                                       
## (Intercept)                                           
## poly(calibrated_cloacal_temp, 2)1                  ***
## poly(calibrated_cloacal_temp, 2)2                  ***
## treatmentCooling                                   ***
## treatmentHeating                                   ***
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling ***
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling ***
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating ***
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.159                                                        
## ply(c__,2)2  0.164 -0.914                                                 
## tretmntClng -0.336  0.260    -0.268                                       
## tretmntHtng -0.324  0.244    -0.251     0.505                             
## pl(__,2)1:C  0.159 -0.999     0.914    -0.259 -0.244                      
## pl(__,2)2:C -0.164  0.914    -1.000     0.268  0.251 -0.913               
## pl(__,2)1:H  0.159 -0.999     0.914    -0.259 -0.244  0.999     -0.913    
## pl(__,2)2:H -0.164  0.914    -1.000     0.268  0.251 -0.914      0.999    
##             p(__,2)1:H
## ply(c__,2)1           
## ply(c__,2)2           
## tretmntClng           
## tretmntHtng           
## pl(__,2)1:C           
## pl(__,2)2:C           
## pl(__,2)1:H           
## pl(__,2)2:H -0.914
#anova(CEWL_LMM_besta, type = "3", ddf = "Kenward-Roger")

CEWL_LMM_bestb <- lmerTest::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             poly(calibrated_dorsal_temp, 2) * treatment +
                             (1|date/lizard_ID))
summary(CEWL_LMM_bestb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 83028.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7628 -0.6345 -0.0833  0.5160  8.3769 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept)  9.996   3.162   
##  date           (Intercept) 21.953   4.685   
##  Residual                    1.554   1.246   
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                                    Estimate Std. Error
## (Intercept)                                           3.248      2.336
## poly(calibrated_dorsal_temp, 2)1                   1534.609     67.430
## poly(calibrated_dorsal_temp, 2)2                  -2711.185     75.106
## treatmentCooling                                      3.529      1.386
## treatmentHeating                                     14.341      1.468
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling -1350.765     67.457
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling  2852.577     75.134
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating -1365.961     67.469
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating  2666.790     75.141
##                                                          df t value Pr(>|t|)
## (Intercept)                                           5.353   1.390   0.2194
## poly(calibrated_dorsal_temp, 2)1                  25237.527  22.759  < 2e-16
## poly(calibrated_dorsal_temp, 2)2                  25236.097 -36.098  < 2e-16
## treatmentCooling                                     30.609   2.546   0.0162
## treatmentHeating                                     30.174   9.766 7.46e-11
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling 25237.523 -20.024  < 2e-16
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling 25236.096  37.967  < 2e-16
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating 25237.517 -20.246  < 2e-16
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating 25236.103  35.490  < 2e-16
##                                                      
## (Intercept)                                          
## poly(calibrated_dorsal_temp, 2)1                  ***
## poly(calibrated_dorsal_temp, 2)2                  ***
## treatmentCooling                                  *  
## treatmentHeating                                  ***
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling ***
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling ***
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating ***
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.162                                                        
## ply(c__,2)2  0.163 -0.946                                                 
## tretmntClng -0.327  0.273    -0.276                                       
## tretmntHtng -0.315  0.257    -0.260     0.507                             
## pl(__,2)1:C  0.162 -1.000     0.946    -0.273 -0.257                      
## pl(__,2)2:C -0.163  0.946    -1.000     0.276  0.259 -0.946               
## pl(__,2)1:H  0.162 -0.999     0.946    -0.273 -0.257  0.999     -0.945    
## pl(__,2)2:H -0.163  0.946    -1.000     0.276  0.259 -0.945      0.999    
##             p(__,2)1:H
## ply(c__,2)1           
## ply(c__,2)2           
## tretmntClng           
## tretmntHtng           
## pl(__,2)1:C           
## pl(__,2)2:C           
## pl(__,2)1:H           
## pl(__,2)2:H -0.945
#anova(CEWL_LMM_bestb, type = "3", ddf = "Kenward-Roger")

# double check sample sizes
temp_time_series %>%
  group_by(lizard_ID, treatment) %>%
  summarise(n()) %>%
  group_by(treatment) %>%
  summarise(n())
## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.
## # A tibble: 3 × 2
##   treatment `n()`
##   <fct>     <int>
## 1 Control      11
## 2 Cooling      12
## 3 Heating      10

These are the best models for CEWL ~ temperature.

broom.mixed::tidy(CEWL_LMM_besta) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/CEWL_best_model_clotemp.csv")

broom.mixed::tidy(CEWL_LMM_bestb) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/CEWL_best_model_dorstemp.csv")

CEWL (raw) ~ time LMM

NOTE: Start times vary from 97-170 seconds after starting time series.

Make a polynomial model first:

rates_LMM_01 <- lme4::lmer(data = CEWL_time_series,
                          CEWL_g_m2h ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
summary(rates_LMM_01)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##    Data: CEWL_time_series
## 
## REML criterion at convergence: 96853.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3325 -0.5700 -0.0438  0.5178  8.0025 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept)  9.976   3.158   
##  date           (Intercept) 17.817   4.221   
##  Residual                    1.754   1.324   
## Number of obs: 28405, groups:  lizard_ID:date, 37; date, 5
## 
## Fixed effects:
##                                     Estimate Std. Error t value
## (Intercept)                           14.957      2.099   7.126
## poly(time_min, 2)1                  -137.974      2.309 -59.765
## poly(time_min, 2)2                    51.100      2.322  22.008
## treatmentCooling                      -8.473      1.297  -6.533
## treatmentHeating                       3.416      1.280   2.669
## poly(time_min, 2)1:treatmentCooling  -88.851      3.376 -26.319
## poly(time_min, 2)2:treatmentCooling   83.367      3.337  24.980
## poly(time_min, 2)1:treatmentHeating  301.172      3.201  94.092
## poly(time_min, 2)2:treatmentHeating -108.483      3.201 -33.893
## 
## Correlation of Fixed Effects:
##             (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1  0.000                                                    
## ply(tm_,2)2  0.000 -0.067                                             
## tretmntClng -0.305  0.001    0.000                                    
## tretmntHtng -0.317  0.001    0.000    0.499                           
## ply(_,2)1:C  0.000 -0.684    0.046    0.000  0.000                    
## ply(_,2)2:C  0.000  0.047   -0.696    0.000  0.000  0.008             
## ply(_,2)1:H  0.000 -0.721    0.049    0.000  0.000  0.493    -0.034   
## ply(_,2)2:H  0.000  0.049   -0.725    0.000  0.000 -0.033     0.505   
##             p(_,2)1:H
## ply(tm_,2)1          
## ply(tm_,2)2          
## tretmntClng          
## tretmntHtng          
## ply(_,2)1:C          
## ply(_,2)2:C          
## ply(_,2)1:H          
## ply(_,2)2:H -0.029
vif(rates_LMM_01)
##                                 GVIF Df GVIF^(1/(2*Df))
## poly(time_min, 2)           9.050897  2        1.734494
## treatment                   1.000002  2        1.000000
## poly(time_min, 2):treatment 9.050904  4        1.317002
drop1(rates_LMM_01)
## Single term deletions
## 
## Model:
## CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##                             npar    AIC
## <none>                            96906
## poly(time_min, 2):treatment    4 111987

Also make a linear model version:

rates_LMM_02 <- lme4::lmer(data = CEWL_time_series,
                          CEWL_g_m2h ~ 
                             time_min * treatment +
                             (1|date/lizard_ID))
summary(rates_LMM_02)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ time_min * treatment + (1 | date/lizard_ID)
##    Data: CEWL_time_series
## 
## REML criterion at convergence: 100915.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6381 -0.5611 -0.0446  0.5129  8.5948 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 10.08    3.175   
##  date           (Intercept) 17.74    4.212   
##  Residual                    2.02    1.421   
## Number of obs: 28405, groups:  lizard_ID:date, 37; date, 5
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)               16.620731   2.097842   7.923
## time_min                  -0.205350   0.003773 -54.431
## treatmentCooling          -7.212880   1.304581  -5.529
## treatmentHeating          -0.277100   1.287525  -0.215
## time_min:treatmentCooling -0.156727   0.005515 -28.418
## time_min:treatmentHeating  0.455474   0.005236  86.986
## 
## Correlation of Fixed Effects:
##             (Intr) tim_mn trtmnC trtmnH tm_m:C
## time_min    -0.015                            
## tretmntClng -0.307  0.024                     
## tretmntHtng -0.319  0.024  0.500              
## tm_mn:trtmC  0.010 -0.684 -0.034 -0.017       
## tm_mn:trtmH  0.011 -0.720 -0.017 -0.033  0.493
anova(rates_LMM_02, rates_LMM_01)
## refitting model(s) with ML (instead of REML)
## Data: CEWL_time_series
## Models:
## rates_LMM_02: CEWL_g_m2h ~ time_min * treatment + (1 | date/lizard_ID)
## rates_LMM_01: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC logLik deviance  Chisq Df Pr(>Chisq)    
## rates_LMM_02    9 100913 100987 -50447   100895                         
## rates_LMM_01   12  96906  97005 -48441    96882 4012.3  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Polynomial is much better than linear.

Also compare to log version:

rates_LMM_03 <- lme4::lmer(data = CEWL_time_series,
                          CEWL_g_m2h ~ 
                             log(time_min) * treatment +
                             (1|date/lizard_ID))
summary(rates_LMM_03)
## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ log(time_min) * treatment + (1 | date/lizard_ID)
##    Data: CEWL_time_series
## 
## REML criterion at convergence: 97113.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1092 -0.5865 -0.0434  0.5019  7.9036 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 10.057   3.171   
##  date           (Intercept) 17.764   4.215   
##  Residual                    1.767   1.329   
## Number of obs: 28405, groups:  lizard_ID:date, 37; date, 5
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    17.83671    2.09858   8.499
## log(time_min)                  -1.48000    0.02353 -62.902
## treatmentCooling               -6.24930    1.30381  -4.793
## treatmentHeating               -2.77380    1.28675  -2.156
## log(time_min):treatmentCooling -1.15769    0.03356 -34.500
## log(time_min):treatmentHeating  3.18416    0.03220  98.902
## 
## Correlation of Fixed Effects:
##             (Intr) lg(t_) trtmnC trtmnH l(_):C
## log(tim_mn) -0.022                            
## tretmntClng -0.307  0.036                     
## tretmntHtng -0.319  0.036  0.500              
## lg(tm_mn):C  0.016 -0.701 -0.050 -0.025       
## lg(tm_mn):H  0.016 -0.731 -0.026 -0.049  0.512
anova(rates_LMM_02, rates_LMM_03)
## refitting model(s) with ML (instead of REML)
## Data: CEWL_time_series
## Models:
## rates_LMM_02: CEWL_g_m2h ~ time_min * treatment + (1 | date/lizard_ID)
## rates_LMM_03: CEWL_g_m2h ~ log(time_min) * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC logLik deviance  Chisq Df Pr(>Chisq)
## rates_LMM_02    9 100913 100987 -50447   100895                     
## rates_LMM_03    9  97121  97196 -48552    97103 3791.4  0
anova(rates_LMM_01, rates_LMM_03)
## refitting model(s) with ML (instead of REML)
## Data: CEWL_time_series
## Models:
## rates_LMM_03: CEWL_g_m2h ~ log(time_min) * treatment + (1 | date/lizard_ID)
## rates_LMM_01: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## rates_LMM_03    9 97121 97196 -48552    97103                         
## rates_LMM_01   12 96906 97005 -48441    96882 220.96  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The polynomial model is best (better than both log and lm options).

Check linear regression assumptions:

plot(rates_LMM_01)

Again, there is a really weird spike, but otherwise there’s no fanning or pattern.

re-run with p-values:

rates_LMM_best <- lmerTest::lmer(data = CEWL_time_series,
                          CEWL_g_m2h ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
summary(rates_LMM_best)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##    Data: CEWL_time_series
## 
## REML criterion at convergence: 96853.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3325 -0.5700 -0.0438  0.5178  8.0025 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept)  9.976   3.158   
##  date           (Intercept) 17.817   4.221   
##  Residual                    1.754   1.324   
## Number of obs: 28405, groups:  lizard_ID:date, 37; date, 5
## 
## Fixed effects:
##                                      Estimate Std. Error        df t value
## (Intercept)                            14.957      2.099     5.189   7.126
## poly(time_min, 2)1                   -137.974      2.309 28362.024 -59.765
## poly(time_min, 2)2                     51.100      2.322 28362.016  22.008
## treatmentCooling                       -8.473      1.297    30.005  -6.533
## treatmentHeating                        3.416      1.280    30.061   2.669
## poly(time_min, 2)1:treatmentCooling   -88.851      3.376 28362.571 -26.319
## poly(time_min, 2)2:treatmentCooling    83.367      3.337 28362.083  24.980
## poly(time_min, 2)1:treatmentHeating   301.172      3.201 28362.093  94.092
## poly(time_min, 2)2:treatmentHeating  -108.483      3.201 28362.048 -33.893
##                                     Pr(>|t|)    
## (Intercept)                         0.000723 ***
## poly(time_min, 2)1                   < 2e-16 ***
## poly(time_min, 2)2                   < 2e-16 ***
## treatmentCooling                    3.17e-07 ***
## treatmentHeating                    0.012151 *  
## poly(time_min, 2)1:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating  < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1  0.000                                                    
## ply(tm_,2)2  0.000 -0.067                                             
## tretmntClng -0.305  0.001    0.000                                    
## tretmntHtng -0.317  0.001    0.000    0.499                           
## ply(_,2)1:C  0.000 -0.684    0.046    0.000  0.000                    
## ply(_,2)2:C  0.000  0.047   -0.696    0.000  0.000  0.008             
## ply(_,2)1:H  0.000 -0.721    0.049    0.000  0.000  0.493    -0.034   
## ply(_,2)2:H  0.000  0.049   -0.725    0.000  0.000 -0.033     0.505   
##             p(_,2)1:H
## ply(tm_,2)1          
## ply(tm_,2)2          
## tretmntClng          
## tretmntHtng          
## ply(_,2)1:C          
## ply(_,2)2:C          
## ply(_,2)1:H          
## ply(_,2)2:H -0.029
#anova(rates_LMM_best, type = "3", ddf = "Kenward-Roger")

# double check sample sizes
CEWL_time_series %>%
  group_by(lizard_ID, treatment) %>%
  summarise(n()) %>%
  group_by(treatment) %>%
  summarise(n())
## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.
## # A tibble: 3 × 2
##   treatment `n()`
##   <fct>     <int>
## 1 Control      12
## 2 Cooling      12
## 3 Heating      13

Treatment explained some variation, and the effect of time explained most of the variation (F(,) = , p <)

broom.mixed::tidy(rates_LMM_best) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/CEWL_best_model_time.csv")

Body Temperature ~ Time

linear models:

body_temp_LMM_01a <- lme4::lmer(data = temp_time_series,
                          calibrated_cloacal_temp ~ 
                             time_min * treatment +
                             (1|date/lizard_ID))
summary(body_temp_LMM_01a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: calibrated_cloacal_temp ~ time_min * treatment + (1 | date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 63649.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9785 -0.4675 -0.0268  0.3484  6.5452 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 3.0775   1.7543  
##  date           (Intercept) 0.1705   0.4129  
##  Residual                   0.7198   0.8484  
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                             Estimate Std. Error  t value
## (Intercept)                26.728519   0.561393   47.611
## time_min                    0.084958   0.002356   36.058
## treatmentCooling            6.075852   0.734617    8.271
## treatmentHeating          -10.116980   0.772484  -13.097
## time_min:treatmentCooling  -1.255658   0.003364 -373.236
## time_min:treatmentHeating   1.161071   0.003411  340.363
## 
## Correlation of Fixed Effects:
##             (Intr) tim_mn trtmnC trtmnH tm_m:C
## time_min    -0.035                            
## tretmntClng -0.680  0.027                     
## tretmntHtng -0.651  0.025  0.491              
## tm_mn:trtmC  0.024 -0.700 -0.037 -0.018       
## tm_mn:trtmH  0.024 -0.691 -0.018 -0.036  0.484
vif(body_temp_LMM_01a)
##                        GVIF Df GVIF^(1/(2*Df))
## time_min           2.875009  1        1.695585
## treatment          1.002620  2        1.000654
## time_min:treatment 2.879911  2        1.302701
body_temp_LMM_01b <- lme4::lmer(data = temp_time_series,
                          calibrated_dorsal_temp ~ 
                             time_min * treatment +
                             (1|date/lizard_ID))
summary(body_temp_LMM_01b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: calibrated_dorsal_temp ~ time_min * treatment + (1 | date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 77541.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8874 -0.3920 -0.0119  0.2875  4.8825 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 1.5354   1.2391  
##  date           (Intercept) 0.3511   0.5925  
##  Residual                   1.2497   1.1179  
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)               27.094037   0.459926   58.91
## time_min                   0.077778   0.003105   25.05
## treatmentCooling           6.418579   0.521574   12.31
## treatmentHeating          -7.961526   0.551634  -14.43
## time_min:treatmentCooling -1.142811   0.004433 -257.81
## time_min:treatmentHeating  0.786695   0.004495  175.02
## 
## Correlation of Fixed Effects:
##             (Intr) tim_mn trtmnC trtmnH tm_m:C
## time_min    -0.056                            
## tretmntClng -0.586  0.050                     
## tretmntHtng -0.562  0.047  0.480              
## tm_mn:trtmC  0.039 -0.700 -0.068 -0.033       
## tm_mn:trtmH  0.039 -0.691 -0.034 -0.066  0.484
vif(body_temp_LMM_01b)
##                        GVIF Df GVIF^(1/(2*Df))
## time_min           2.875083  1        1.695607
## treatment          1.008948  2        1.002230
## time_min:treatment 2.891679  2        1.304030

polynomial models:

body_temp_LMM_02a <- lme4::lmer(data = temp_time_series,
                          calibrated_cloacal_temp ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
summary(body_temp_LMM_02a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: calibrated_cloacal_temp ~ poly(time_min, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 41254.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8989 -0.4649 -0.0011  0.3942  6.3225 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 3.0570   1.7484  
##  date           (Intercept) 0.1856   0.4308  
##  Residual                   0.2967   0.5447  
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                      Estimate Std. Error  t value
## (Intercept)                           27.4181     0.5621   48.780
## poly(time_min, 2)1                    52.0361     0.9376   55.500
## poly(time_min, 2)2                     6.0643     0.9435    6.428
## treatmentCooling                      -4.0686     0.7317   -5.560
## treatmentHeating                      -0.7262     0.7697   -0.944
## poly(time_min, 2)1:treatmentCooling -764.1781     1.3382 -571.054
## poly(time_min, 2)2:treatmentCooling  161.4913     1.3255  121.838
## poly(time_min, 2)1:treatmentHeating  716.1235     1.3557  528.233
## poly(time_min, 2)2:treatmentHeating  -63.6826     1.3531  -47.063
## 
## Correlation of Fixed Effects:
##             (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001                                                    
## ply(tm_,2)2  0.000 -0.077                                             
## tretmntClng -0.676  0.000    0.000                                    
## tretmntHtng -0.647  0.000    0.000    0.490                           
## ply(_,2)1:C  0.000 -0.701    0.054    0.000  0.000                    
## ply(_,2)2:C  0.000  0.055   -0.712    0.000  0.000 -0.005             
## ply(_,2)1:H  0.000 -0.692    0.053    0.000  0.000  0.485    -0.038   
## ply(_,2)2:H  0.000  0.054   -0.697    0.000  0.000 -0.038     0.496   
##             p(_,2)1:H
## ply(tm_,2)1          
## ply(tm_,2)2          
## tretmntClng          
## tretmntHtng          
## ply(_,2)1:C          
## ply(_,2)2:C          
## ply(_,2)1:H          
## ply(_,2)2:H -0.021
car::vif(body_temp_LMM_02a)
##                                 GVIF Df GVIF^(1/(2*Df))
## poly(time_min, 2)           8.580314  2        1.711496
## treatment                   1.000001  2        1.000000
## poly(time_min, 2):treatment 8.580316  4        1.308241
anova(body_temp_LMM_02a, body_temp_LMM_01a)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## body_temp_LMM_01a: calibrated_cloacal_temp ~ time_min * treatment + (1 | date/lizard_ID)
## body_temp_LMM_02a: calibrated_cloacal_temp ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##                   npar   AIC   BIC logLik deviance Chisq Df Pr(>Chisq)    
## body_temp_LMM_01a    9 63638 63712 -31810    63620                        
## body_temp_LMM_02a   12 41291 41388 -20633    41267 22354  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
body_temp_LMM_02b <- lme4::lmer(data = temp_time_series,
                          calibrated_dorsal_temp ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
summary(body_temp_LMM_02b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: calibrated_dorsal_temp ~ poly(time_min, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 71392.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8879 -0.3478 -0.0225  0.3300  6.0970 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 1.5384   1.2403  
##  date           (Intercept) 0.3563   0.5969  
##  Residual                   0.9812   0.9906  
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                      Estimate Std. Error  t value
## (Intercept)                           27.7236     0.4606   60.191
## poly(time_min, 2)1                    48.0115     1.7051   28.158
## poly(time_min, 2)2                     0.6860     1.7158    0.400
## treatmentCooling                      -2.8157     0.5208   -5.406
## treatmentHeating                      -1.5969     0.5509   -2.899
## poly(time_min, 2)1:treatmentCooling -697.0769     2.4336 -286.444
## poly(time_min, 2)2:treatmentCooling  133.9832     2.4104   55.585
## poly(time_min, 2)1:treatmentHeating  484.8482     2.4654  196.660
## poly(time_min, 2)2:treatmentHeating  -42.9964     2.4608  -17.473
## 
## Correlation of Fixed Effects:
##             (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001                                                    
## ply(tm_,2)2  0.000 -0.077                                             
## tretmntClng -0.584  0.001    0.000                                    
## tretmntHtng -0.560  0.001    0.000    0.480                           
## ply(_,2)1:C  0.001 -0.701    0.054    0.000 -0.001                    
## ply(_,2)2:C  0.000  0.055   -0.712    0.000  0.000 -0.005             
## ply(_,2)1:H  0.001 -0.692    0.053   -0.001  0.000  0.485    -0.038   
## ply(_,2)2:H  0.000  0.054   -0.697    0.000  0.000 -0.038     0.496   
##             p(_,2)1:H
## ply(tm_,2)1          
## ply(tm_,2)2          
## tretmntClng          
## tretmntHtng          
## ply(_,2)1:C          
## ply(_,2)2:C          
## ply(_,2)1:H          
## ply(_,2)2:H -0.021
car::vif(body_temp_LMM_02b)
##                                 GVIF Df GVIF^(1/(2*Df))
## poly(time_min, 2)           8.580549  2        1.711507
## treatment                   1.000007  2        1.000002
## poly(time_min, 2):treatment 8.580567  4        1.308246
anova(body_temp_LMM_02b, body_temp_LMM_01b)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## body_temp_LMM_01b: calibrated_dorsal_temp ~ time_min * treatment + (1 | date/lizard_ID)
## body_temp_LMM_02b: calibrated_dorsal_temp ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##                   npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## body_temp_LMM_01b    9 77531 77604 -38756    77513                         
## body_temp_LMM_02b   12 71434 71532 -35705    71410 6102.6  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Again, the polynomial models are MUCH better.

body_temp_LMM_besta <- lmerTest::lmer(data = temp_time_series,
                          calibrated_cloacal_temp ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
summary(body_temp_LMM_besta)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: calibrated_cloacal_temp ~ poly(time_min, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 41254.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8989 -0.4649 -0.0011  0.3942  6.3225 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 3.0570   1.7484  
##  date           (Intercept) 0.1856   0.4308  
##  Residual                   0.2967   0.5447  
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                       Estimate Std. Error         df  t value
## (Intercept)                            27.4181     0.5621    16.5563   48.780
## poly(time_min, 2)1                     52.0361     0.9376 25212.0138   55.500
## poly(time_min, 2)2                      6.0643     0.9435 25212.0089    6.428
## treatmentCooling                       -4.0686     0.7317    26.4065   -5.560
## treatmentHeating                       -0.7262     0.7697    27.6621   -0.944
## poly(time_min, 2)1:treatmentCooling  -764.1781     1.3382 25212.2800 -571.054
## poly(time_min, 2)2:treatmentCooling   161.4913     1.3255 25212.0407  121.838
## poly(time_min, 2)1:treatmentHeating   716.1235     1.3557 25212.0589  528.233
## poly(time_min, 2)2:treatmentHeating   -63.6826     1.3531 25212.0326  -47.063
##                                     Pr(>|t|)    
## (Intercept)                          < 2e-16 ***
## poly(time_min, 2)1                   < 2e-16 ***
## poly(time_min, 2)2                  1.32e-10 ***
## treatmentCooling                    7.32e-06 ***
## treatmentHeating                       0.354    
## poly(time_min, 2)1:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating  < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001                                                    
## ply(tm_,2)2  0.000 -0.077                                             
## tretmntClng -0.676  0.000    0.000                                    
## tretmntHtng -0.647  0.000    0.000    0.490                           
## ply(_,2)1:C  0.000 -0.701    0.054    0.000  0.000                    
## ply(_,2)2:C  0.000  0.055   -0.712    0.000  0.000 -0.005             
## ply(_,2)1:H  0.000 -0.692    0.053    0.000  0.000  0.485    -0.038   
## ply(_,2)2:H  0.000  0.054   -0.697    0.000  0.000 -0.038     0.496   
##             p(_,2)1:H
## ply(tm_,2)1          
## ply(tm_,2)2          
## tretmntClng          
## tretmntHtng          
## ply(_,2)1:C          
## ply(_,2)2:C          
## ply(_,2)1:H          
## ply(_,2)2:H -0.021
write.csv(broom.mixed::tidy(body_temp_LMM_besta), 
          "./Results_Statistics/temp_time_best_model_clotemp.csv")

body_temp_LMM_bestb <- lmerTest::lmer(data = temp_time_series,
                          calibrated_dorsal_temp ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
summary(body_temp_LMM_bestb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: calibrated_dorsal_temp ~ poly(time_min, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 71392.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8879 -0.3478 -0.0225  0.3300  6.0970 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 1.5384   1.2403  
##  date           (Intercept) 0.3563   0.5969  
##  Residual                   0.9812   0.9906  
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                       Estimate Std. Error         df  t value
## (Intercept)                            27.7236     0.4606    11.4180   60.191
## poly(time_min, 2)1                     48.0115     1.7051 25212.0897   28.158
## poly(time_min, 2)2                      0.6860     1.7158 25212.0585    0.400
## treatmentCooling                       -2.8157     0.5208    26.3868   -5.406
## treatmentHeating                       -1.5969     0.5509    27.0887   -2.899
## poly(time_min, 2)1:treatmentCooling  -697.0769     2.4336 25213.7698 -286.444
## poly(time_min, 2)2:treatmentCooling   133.9832     2.4104 25212.2645   55.585
## poly(time_min, 2)1:treatmentHeating   484.8482     2.4654 25212.3649  196.660
## poly(time_min, 2)2:treatmentHeating   -42.9964     2.4608 25212.2024  -17.473
##                                     Pr(>|t|)    
## (Intercept)                         1.19e-15 ***
## poly(time_min, 2)1                   < 2e-16 ***
## poly(time_min, 2)2                   0.68930    
## treatmentCooling                    1.10e-05 ***
## treatmentHeating                     0.00734 ** 
## poly(time_min, 2)1:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating  < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001                                                    
## ply(tm_,2)2  0.000 -0.077                                             
## tretmntClng -0.584  0.001    0.000                                    
## tretmntHtng -0.560  0.001    0.000    0.480                           
## ply(_,2)1:C  0.001 -0.701    0.054    0.000 -0.001                    
## ply(_,2)2:C  0.000  0.055   -0.712    0.000  0.000 -0.005             
## ply(_,2)1:H  0.001 -0.692    0.053   -0.001  0.000  0.485    -0.038   
## ply(_,2)2:H  0.000  0.054   -0.697    0.000  0.000 -0.038     0.496   
##             p(_,2)1:H
## ply(tm_,2)1          
## ply(tm_,2)2          
## tretmntClng          
## tretmntHtng          
## ply(_,2)1:C          
## ply(_,2)2:C          
## ply(_,2)1:H          
## ply(_,2)2:H -0.021
write.csv(broom.mixed::tidy(body_temp_LMM_bestb), 
          "./Results_Statistics/temp_time_best_model_dorstemp.csv")

Rate of Change

We want to know why each lizard had a different rate of change of CEWL during the 15-minute temperature treatment and measurement period.

First, calculate rate of change for each lizard:

change_LM <- lm(data = CEWL_time_series,
                          CEWL_g_m2h ~ 
                             time_min*lizard_ID)
summary(change_LM)
## 
## Call:
## lm(formula = CEWL_g_m2h ~ time_min * lizard_ID, data = CEWL_time_series)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8164 -0.5799 -0.0359  0.4370  9.2347 
## 
## Coefficients:
##                         Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)            26.849726   0.079132  339.303  < 2e-16 ***
## time_min                0.122080   0.008697   14.037  < 2e-16 ***
## lizard_ID402           -0.391770   0.111914   -3.501 0.000465 ***
## lizard_ID404           -4.828469   0.110317  -43.769  < 2e-16 ***
## lizard_ID406          -15.285763   0.118305 -129.207  < 2e-16 ***
## lizard_ID407           -0.098935   0.109533   -0.903 0.366404    
## lizard_ID408           -1.653538   0.122070  -13.546  < 2e-16 ***
## lizard_ID409           -6.242781   0.113576  -54.966  < 2e-16 ***
## lizard_ID410          -24.535008   0.111329 -220.384  < 2e-16 ***
## lizard_ID411          -17.712482   0.115498 -153.358  < 2e-16 ***
## lizard_ID412          -16.092335   0.110917 -145.084  < 2e-16 ***
## lizard_ID413           -8.859382   0.113889  -77.789  < 2e-16 ***
## lizard_ID414          -21.472142   0.126188 -170.159  < 2e-16 ***
## lizard_ID415          -17.170442   0.115574 -148.567  < 2e-16 ***
## lizard_ID416           -6.509572   0.111739  -58.257  < 2e-16 ***
## lizard_ID417          -22.062135   0.111838 -197.269  < 2e-16 ***
## lizard_ID418          -18.675743   0.111729 -167.152  < 2e-16 ***
## lizard_ID419          -16.585398   0.112004 -148.079  < 2e-16 ***
## lizard_ID420          -21.060422   0.112011 -188.022  < 2e-16 ***
## lizard_ID421          -20.484787   0.112438 -182.188  < 2e-16 ***
## lizard_ID422          -14.092361   0.114180 -123.422  < 2e-16 ***
## lizard_ID423          -21.460844   0.111367 -192.703  < 2e-16 ***
## lizard_ID424          -19.003214   0.111240 -170.831  < 2e-16 ***
## lizard_ID425          -16.936058   0.111112 -152.424  < 2e-16 ***
## lizard_ID426          -19.507882   0.148579 -131.296  < 2e-16 ***
## lizard_ID427          -14.604353   0.113202 -129.012  < 2e-16 ***
## lizard_ID428           -3.577002   0.120484  -29.689  < 2e-16 ***
## lizard_ID429          -12.050071   0.111385 -108.183  < 2e-16 ***
## lizard_ID430          -13.697643   0.112326 -121.945  < 2e-16 ***
## lizard_ID431          -18.738456   0.111911 -167.441  < 2e-16 ***
## lizard_ID432          -19.967832   0.112583 -177.360  < 2e-16 ***
## lizard_ID433           -9.059766   0.112296  -80.677  < 2e-16 ***
## lizard_ID434           -5.041610   0.112862  -44.671  < 2e-16 ***
## lizard_ID435          -12.863750   0.112050 -114.804  < 2e-16 ***
## lizard_ID436           -5.219401   0.123891  -42.129  < 2e-16 ***
## lizard_ID437           -8.399519   0.112705  -74.527  < 2e-16 ***
## lizard_ID438           -6.781302   0.113911  -59.532  < 2e-16 ***
## lizard_ID439          -11.307001   0.111506 -101.403  < 2e-16 ***
## time_min:lizard_ID402  -0.385888   0.012300  -31.373  < 2e-16 ***
## time_min:lizard_ID404  -0.494876   0.012165  -40.679  < 2e-16 ***
## time_min:lizard_ID406  -0.793039   0.013514  -58.683  < 2e-16 ***
## time_min:lizard_ID407   0.100523   0.012103    8.306  < 2e-16 ***
## time_min:lizard_ID408  -0.561813   0.013056  -43.031  < 2e-16 ***
## time_min:lizard_ID409  -0.622057   0.012433  -50.033  < 2e-16 ***
## time_min:lizard_ID410  -0.020582   0.012247   -1.681 0.092869 .  
## time_min:lizard_ID411  -0.463220   0.012420  -37.295  < 2e-16 ***
## time_min:lizard_ID412   0.237357   0.012215   19.432  < 2e-16 ***
## time_min:lizard_ID413  -0.341593   0.012458  -27.419  < 2e-16 ***
## time_min:lizard_ID414  -0.209479   0.017039  -12.294  < 2e-16 ***
## time_min:lizard_ID415  -0.205661   0.012599  -16.324  < 2e-16 ***
## time_min:lizard_ID416  -0.377975   0.012282  -30.775  < 2e-16 ***
## time_min:lizard_ID417   0.291636   0.012292   23.725  < 2e-16 ***
## time_min:lizard_ID418   0.211037   0.012282   17.182  < 2e-16 ***
## time_min:lizard_ID419   0.118183   0.012305    9.604  < 2e-16 ***
## time_min:lizard_ID420  -0.340719   0.012310  -27.679  < 2e-16 ***
## time_min:lizard_ID421  -0.407892   0.012343  -33.047  < 2e-16 ***
## time_min:lizard_ID422  -0.293826   0.012488  -23.529  < 2e-16 ***
## time_min:lizard_ID423   0.567593   0.012255   46.313  < 2e-16 ***
## time_min:lizard_ID424   0.573931   0.012241   46.888  < 2e-16 ***
## time_min:lizard_ID425  -0.183705   0.012230  -15.021  < 2e-16 ***
## time_min:lizard_ID426  -0.820460   0.025602  -32.047  < 2e-16 ***
## time_min:lizard_ID427  -0.013238   0.012403   -1.067 0.285850    
## time_min:lizard_ID428  -0.157217   0.012993  -12.100  < 2e-16 ***
## time_min:lizard_ID429  -0.061577   0.012255   -5.025 5.07e-07 ***
## time_min:lizard_ID430  -0.924418   0.012332  -74.963  < 2e-16 ***
## time_min:lizard_ID431  -0.440182   0.012298  -35.793  < 2e-16 ***
## time_min:lizard_ID432  -0.284201   0.012358  -22.998  < 2e-16 ***
## time_min:lizard_ID433  -0.639850   0.012329  -51.897  < 2e-16 ***
## time_min:lizard_ID434  -0.513333   0.012377  -41.473  < 2e-16 ***
## time_min:lizard_ID435  -0.679765   0.012310  -55.219  < 2e-16 ***
## time_min:lizard_ID436  -0.037103   0.015879   -2.337 0.019469 *  
## time_min:lizard_ID437  -0.574387   0.012365  -46.452  < 2e-16 ***
## time_min:lizard_ID438  -0.282519   0.012464  -22.667  < 2e-16 ***
## time_min:lizard_ID439   0.538159   0.012266   43.874  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9663 on 28331 degrees of freedom
## Multiple R-squared:  0.9816, Adjusted R-squared:  0.9815 
## F-statistic: 2.069e+04 on 73 and 28331 DF,  p-value: < 2.2e-16
# now get the trend/slope for each individual and make into a variable in df
lizard_coeffs <- data.frame(emtrends(change_LM, "lizard_ID", var = "time_min")) %>%
   left_join(read_rds("Data/collection_dat_formatted.RDS"), 
             by = c("lizard_ID" = "individual_ID")) %>%
   dplyr::select(lizard_ID, rate_change = time_min.trend, 
                 SE, df, lower.CL, upper.CL,
                 date,
                 treatment,
                 mass_g, SVL_mm,
                 chamber_time_elapsed_hr
                 ) %>%
   dplyr::mutate(date_factor = as.factor(date))
summary(lizard_coeffs)
##    lizard_ID   rate_change            SE                 df       
##  401    : 1   Min.   :-0.8023   Min.   :0.008416   Min.   :28331  
##  402    : 1   1st Qu.:-0.3728   1st Qu.:0.008672   1st Qu.:28331  
##  404    : 1   Median :-0.1621   Median :0.008739   Median :28331  
##  406    : 1   Mean   :-0.1074   Mean   :0.009524   Mean   :28331  
##  407    : 1   3rd Qu.: 0.1088   3rd Qu.:0.008920   3rd Qu.:28331  
##  408    : 1   Max.   : 0.6960   Max.   :0.024080   Max.   :28331  
##  (Other):31                                                       
##     lower.CL           upper.CL             date              treatment 
##  Min.   :-0.81947   Min.   :-0.78520   Min.   :2021-10-05   Control:12  
##  1st Qu.:-0.38947   1st Qu.:-0.35612   1st Qu.:2021-10-11   Cooling:12  
##  Median :-0.17933   Median :-0.14491   Median :2021-10-12   Heating:13  
##  Mean   :-0.12610   Mean   :-0.08877   Mean   :2021-10-13               
##  3rd Qu.: 0.09151   3rd Qu.: 0.12617   3rd Qu.:2021-10-18               
##  Max.   : 0.67913   Max.   : 0.71289   Max.   :2021-10-19               
##                                                                         
##      mass_g          SVL_mm      chamber_time_elapsed_hr     date_factor
##  Min.   : 8.70   Min.   :60.00   Min.   :0.9333          2021-10-05:6   
##  1st Qu.:10.00   1st Qu.:64.00   1st Qu.:1.0000          2021-10-11:8   
##  Median :11.30   Median :66.00   Median :1.0333          2021-10-12:7   
##  Mean   :11.42   Mean   :66.24   Mean   :1.1005          2021-10-18:8   
##  3rd Qu.:12.60   3rd Qu.:69.00   3rd Qu.:1.1500          2021-10-19:8   
##  Max.   :15.30   Max.   :75.00   Max.   :1.5333                         
## 
length(unique(CEWL_time_series$lizard_ID))
## [1] 37

NOW we can look at how the RATE of CEWL change (CEWL change per MINUTE) is different among lizards.

change_LMM_1 <- lme4::lmer(data = lizard_coeffs,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                       mass_g + SVL_mm + 
                       chamber_time_elapsed_hr + 
                  # include date as random effect bc sig diff weather
                       (1|date_factor))
summary(change_LMM_1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: rate_change ~ treatment + mass_g + SVL_mm + chamber_time_elapsed_hr +  
##     (1 | date_factor)
##    Data: lizard_coeffs
## 
## REML criterion at convergence: 26.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.97945 -0.60829  0.04974  0.48861  1.63798 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  date_factor (Intercept) 0.01674  0.1294  
##  Residual                0.07171  0.2678  
## Number of obs: 37, groups:  date_factor, 5
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)             -0.90245    1.03318  -0.873
## treatmentCooling        -0.12485    0.14578  -0.856
## treatmentHeating         0.49210    0.11557   4.258
## mass_g                  -0.03169    0.03778  -0.839
## SVL_mm                   0.01927    0.01781   1.082
## chamber_time_elapsed_hr -0.23232    0.42633  -0.545
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnC trtmnH mass_g SVL_mm
## tretmntClng  0.194                            
## tretmntHtng -0.165  0.533                     
## mass_g       0.396  0.060 -0.084              
## SVL_mm      -0.880  0.011  0.227 -0.650       
## chmbr_tm_l_ -0.346 -0.656 -0.271 -0.157 -0.002
car::vif(change_LMM_1)
##                             GVIF Df GVIF^(1/(2*Df))
## treatment               1.934248  2        1.179310
## mass_g                  1.825412  1        1.351078
## SVL_mm                  1.878111  1        1.370442
## chamber_time_elapsed_hr 1.865277  1        1.365751
drop1(change_LMM_1)
## Single term deletions
## 
## Model:
## rate_change ~ treatment + mass_g + SVL_mm + chamber_time_elapsed_hr + 
##     (1 | date_factor)
##                         npar    AIC
## <none>                       21.929
## treatment                  2 42.116
## mass_g                     1 20.523
## SVL_mm                     1 21.118
## chamber_time_elapsed_hr    1 20.291
anova(change_LMM_1)
## Analysis of Variance Table
##                         npar  Sum Sq Mean Sq F value
## treatment                  2 2.56693 1.28346 17.8977
## mass_g                     1 0.00639 0.00639  0.0891
## SVL_mm                     1 0.08371 0.08371  1.1673
## chamber_time_elapsed_hr    1 0.02129 0.02129  0.2969

Don’t need to worry about VIF.

drop mass first based on low SS and resulting AIC:

change_LMM_2 <- lme4::lmer(data = lizard_coeffs,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                       SVL_mm + 
                       chamber_time_elapsed_hr + 
                  # include date as random effect bc sig diff weather
                       (1|date_factor))
summary(change_LMM_2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: rate_change ~ treatment + SVL_mm + chamber_time_elapsed_hr +  
##     (1 | date_factor)
##    Data: lizard_coeffs
## 
## REML criterion at convergence: 22.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8071 -0.6048  0.1496  0.5868  1.7654 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  date_factor (Intercept) 0.01111  0.1054  
##  Residual                0.07302  0.2702  
## Number of obs: 37, groups:  date_factor, 5
## 
## Fixed effects:
##                          Estimate Std. Error t value
## (Intercept)             -0.564781   0.948384  -0.596
## treatmentCooling        -0.114887   0.146720  -0.783
## treatmentHeating         0.483926   0.116036   4.170
## SVL_mm                   0.009684   0.013537   0.715
## chamber_time_elapsed_hr -0.290802   0.424397  -0.685
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnC trtmnH SVL_mm
## tretmntClng  0.183                     
## tretmntHtng -0.139  0.544              
## SVL_mm      -0.891  0.070  0.224       
## chmbr_tm_l_ -0.315 -0.656 -0.291 -0.139
drop1(change_LMM_2)
## Single term deletions
## 
## Model:
## rate_change ~ treatment + SVL_mm + chamber_time_elapsed_hr + 
##     (1 | date_factor)
##                         npar    AIC
## <none>                       20.523
## treatment                  2 40.292
## SVL_mm                     1 19.118
## chamber_time_elapsed_hr    1 19.059
anova(change_LMM_2)
## Analysis of Variance Table
##                         npar  Sum Sq Mean Sq F value
## treatment                  2 2.56389 1.28195 17.5558
## SVL_mm                     1 0.02860 0.02860  0.3917
## chamber_time_elapsed_hr    1 0.03428 0.03428  0.4695

drop SVL:

change_LMM_3 <- lme4::lmer(data = lizard_coeffs,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                             chamber_time_elapsed_hr + 
                  # include date as random effect bc sig diff weather
                       (1|date_factor))
summary(change_LMM_3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: rate_change ~ treatment + chamber_time_elapsed_hr + (1 | date_factor)
##    Data: lizard_coeffs
## 
## REML criterion at convergence: 16.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8677 -0.6281  0.1374  0.5128  1.7759 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  date_factor (Intercept) 0.01147  0.1071  
##  Residual                0.07171  0.2678  
## Number of obs: 37, groups:  date_factor, 5
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)              0.03945    0.42710   0.092
## treatmentCooling        -0.12250    0.14505  -0.845
## treatmentHeating         0.46534    0.11207   4.152
## chamber_time_elapsed_hr -0.24832    0.41652  -0.596
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnC trtmnH
## tretmntClng  0.543              
## tretmntHtng  0.136  0.543       
## chmbr_tm_l_ -0.977 -0.654 -0.269
drop1(change_LMM_3)
## Single term deletions
## 
## Model:
## rate_change ~ treatment + chamber_time_elapsed_hr + (1 | date_factor)
##                         npar    AIC
## <none>                       19.118
## treatment                  2 38.317
## chamber_time_elapsed_hr    1 17.510
anova(change_LMM_3)
## Analysis of Variance Table
##                         npar  Sum Sq Mean Sq F value
## treatment                  2 2.56421 1.28211 17.8795
## chamber_time_elapsed_hr    1 0.02549 0.02549  0.3554

drop chamber_time_elapsed_hr:

change_LMM_4 <- lme4::lmer(data = lizard_coeffs,
                          rate_change ~ 
                             treatment + 
                       (1|date_factor))
summary(change_LMM_4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: rate_change ~ treatment + (1 | date_factor)
##    Data: lizard_coeffs
## 
## REML criterion at convergence: 16.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.94436 -0.62201 -0.04392  0.39606  1.77788 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  date_factor (Intercept) 0.01179  0.1086  
##  Residual                0.07012  0.2648  
## Number of obs: 37, groups:  date_factor, 5
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       -0.2093     0.0908  -2.305
## treatmentCooling  -0.1793     0.1085  -1.653
## treatmentHeating   0.4475     0.1068   4.191
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnC
## tretmntClng -0.593       
## tretmntHtng -0.611  0.504
drop1(change_LMM_4)
## boundary (singular) fit: see help('isSingular')
## Single term deletions
## 
## Model:
## rate_change ~ treatment + (1 | date_factor)
##           npar    AIC
## <none>         17.510
## treatment    2 38.876
anova(change_LMM_4)
## Analysis of Variance Table
##           npar Sum Sq Mean Sq F value
## treatment    2 2.5645  1.2823  18.286
plot(change_LMM_4)

The simplest model is our best/final model, and assumptions look good.

Save with p-values:

change_LMM_best <- lmerTest::lmer(data = lizard_coeffs,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                       (1|date_factor))
summary(change_LMM_best)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rate_change ~ treatment + (1 | date_factor)
##    Data: lizard_coeffs
## 
## REML criterion at convergence: 16.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.94436 -0.62201 -0.04392  0.39606  1.77788 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  date_factor (Intercept) 0.01179  0.1086  
##  Residual                0.07012  0.2648  
## Number of obs: 37, groups:  date_factor, 5
## 
## Fixed effects:
##                  Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)       -0.2093     0.0908 12.1290  -2.305 0.039599 *  
## treatmentCooling  -0.1793     0.1085 29.9201  -1.653 0.108748    
## treatmentHeating   0.4475     0.1068 30.3121   4.191 0.000222 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnC
## tretmntClng -0.593       
## tretmntHtng -0.611  0.504
broom.mixed::tidy(change_LMM_best) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/change_CEWL_best_model.csv")

anova(change_LMM_best, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##           Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## treatment 2.5102  1.2551     2 30.708  17.898 7.034e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Treatment group explained a significant amount of the variation (F(3,30)=18.29, p < 0.0001).

also get pairwise differences for change ~ tmt, which I will use to make the boxplot later:

emmeans <- data.frame(emmeans(change_LMM_best, "treatment"))
emmeans
##   treatment     emmean         SE       df    lower.CL    upper.CL
## 1   Control -0.2093260 0.09102635 13.09716 -0.40582833 -0.01282375
## 2   Cooling -0.3886322 0.09189944 12.94683 -0.58725181 -0.19001265
## 3   Heating  0.2381211 0.08845166 11.97139  0.04535039  0.43089180
emmean_ps <- data.frame(pairs(emmeans(change_LMM_best, "treatment"))) # p values
emmean_ps
##            contrast   estimate        SE       df   t.ratio      p.value
## 1 Control - Cooling  0.1793062 0.1087886 30.33842  1.648208 2.414730e-01
## 2 Control - Heating -0.4474471 0.1074827 30.69471 -4.162968 6.708347e-04
## 3 Cooling - Heating -0.6267533 0.1082740 31.09322 -5.788584 6.511073e-06

Figures

COLOR

cool <- c(brewer.pal(11, "RdBu")[c(11)])
heat <- c(brewer.pal(11, "RdBu")[c(2)])
control <- c(brewer.pal(9, "Greys")[c(5)])
cntrls <- c(brewer.pal(9, "Greys")[c(3,6,9,4,7,3,5,8,4,6,9)])

cloacal temp ~ time

ggplot(temp_time_series) +
   aes(x = time_sec, 
       y = calibrated_cloacal_temp,
       color = treatment) +
   geom_point(size = 0.01, 
              alpha = 0.2) +
   geom_smooth(method = "loess", 
               se = F, 
               size = 2,
               na.rm = TRUE) +
   theme_classic() +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time (seconds)" , 
        y = "Internal Body Temperature (°C)", 
        color = "Treatment", se = FALSE) -> ctemp_time
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
ctemp_time
## `geom_smooth()` using formula = 'y ~ x'

dorsal temp ~ time

ggplot(temp_time_series) +
   aes(x = time_sec, 
       y = (calibrated_dorsal_temp),
       color = treatment) +
   geom_point(size = 0.01, 
              alpha = 0.2) +
   geom_smooth(method = "loess", 
               se = F, 
               size = 2,
               na.rm = TRUE) +
   theme_classic() +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x= "Time (seconds)" , 
        y= "Surface Body Temperature (°C)", 
        color ="Treatment") -> dtemp_time
dtemp_time
## `geom_smooth()` using formula = 'y ~ x'

CEWL (raw) ~ time

ggplot(CEWL_time_series) +
   aes(x= time_sec, 
       y = CEWL_g_m2h,
       color = treatment) +
   geom_point(size = 0.01, 
              alpha = 0.25) +
   geom_smooth(method = "loess", 
               size = 2.5, 
               se = F) +
   theme_classic() +
   theme(text = element_text(size = 15)) + 
   theme(axis.text.x = element_text(size = 10)) +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time (seconds)" , 
        y = bquote('CEWL (g/'*m^2*'h)'), 
        color = "Treatment", 
        se = FALSE) +
   scale_x_continuous(limits = c(0, 925), 
                      breaks=seq(0, 950, 100)) -> CEWL_time
CEWL_time
## `geom_smooth()` using formula = 'y ~ x'

legend

ggplot(temp_time_series) +
   
   geom_smooth(aes(x = (calibrated_cloacal_temp), 
       y = CEWL_g_m2h,
       color = treatment),
       method = "lm", 
        formula = y ~ poly(x, 2),
               size = 1, 
               se = F) +
  geom_point(aes(x = (calibrated_cloacal_temp), 
       y = CEWL_g_m2h,
       fill = treatment,
       shape = treatment),
       size = 4, 
       alpha = 1) +
  
  scale_x_continuous(limits = c(14, 39),
                     breaks = c(seq(15, 35, 5)),
                     labels = c(seq(15, 35, 5))) +
  theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = c(0.5,0.5)) +
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") +
   scale_fill_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  scale_color_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment")  -> fancy_legend
  
  
fancy_legend

LEGEND <- as_ggplot(get_legend(fancy_legend)) 

# save for using in arranged multi-plots
ggsave(filename = "legend_only.pdf",
       plot = LEGEND,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 40, height = 30)

CEWL (raw) ~ cloacal temperature

ggplot(temp_time_series) +
   geom_point(aes(x = (calibrated_cloacal_temp), 
       y = CEWL_g_m2h,
       shape = treatment,
       color = treatment),
       size = 0.01, 
               alpha = 0.1) +
   geom_smooth(aes(x = (calibrated_cloacal_temp), 
       y = CEWL_g_m2h,
       color = treatment),
       method = "lm", 
        formula = y ~ poly(x, 2),
               size = 1.5, 
               se = F) +
  scale_x_continuous(limits = c(12.5, 40),
                     breaks = c(seq(15, 35, 5)),
                     labels = c(seq(15, 35, 5))) +

  # control group
   #annotate(geom = "point", x = 25.4, y = 9.6, 
    #       size = 5, pch = 21, fill = control) + 
   #annotate(geom = "point", x = 30.14, y = 21.4, 
    #       size = 5, pch = 21, fill = control) + 
  
  # annotations for cooling group
   annotate(geom = "point", x = 14.3, y = 5, 
           size = 4, pch = 25, fill = cool) + 
   annotate("text", x = 12.6, y = 5, label = "T[15]", parse = T,
           size = 5, family = "sans", color = cool) +
   annotate(geom = "point", x = 35.2, y = 11.8, 
           size = 4, pch = 25, fill = cool) + 
   annotate("text", x = 36.7, y = 11.8, label = "T[1]", parse = T,
           size = 5, family = "sans", color = cool) +
  
  # annotations for heating group
   annotate(geom = "point", x = 15.9, y = 14, 
           size = 4, pch = 24, fill = heat) +  
   annotate("text", x = 14.3, y = 14, label = "T[1]", parse = T,
           size = 5, family = "sans", color = heat) +
   annotate(geom = "point", x = 38.3, y = 25.2, 
           size = 4, pch = 24, fill = heat) +
   annotate("text", x = 40, y = 25.2, label = "T[15]", parse = T,
           size = 5, family = "sans", color = heat) +
   
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        #legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
  
   scale_color_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
  
   labs(x = "Internal Body Temperature (°C)" , 
        y = bquote('CEWL (g '*m^-2*' '*h^-1*')'), 
        se = FALSE) -> CEWL_ctemp
CEWL_ctemp

# save
ggsave(filename = "CEWL_int_body_temp.pdf",
       plot = CEWL_ctemp,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 150, height = 100)

CEWL Control

# get control data only
temp_time_series %>%
  dplyr::filter(treatment == "Control") %>%
  dplyr::filter(complete.cases(calibrated_cloacal_temp, CEWL_g_m2h)) %>%
  group_by(lizard_ID) %>%
  summarise(max(CEWL_g_m2h)) %>%
  nrow()
## [1] 11
# plot control data only
temp_time_series %>%
  dplyr::filter(treatment == "Control") %>%
  ggplot() +
   aes(x = calibrated_cloacal_temp, 
       y = CEWL_g_m2h,
       color = lizard_ID) +
  
   geom_point(size = 0.1, 
               alpha = 0.1) +
   geom_smooth(method = "lm", 
               size = 1, 
               se = F) +

   ylim(0,32) +
   xlim(25,31) +
  
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
  
   scale_color_manual(values = cntrls) +
  
   labs(x = "Internal Body Temperature (°C)" , 
        y = bquote('CEWL (g '*m^-2*' '*h^-1*')'), 
        color = "Treatment", 
        se = FALSE) -> CEWL_control
CEWL_control
## `geom_smooth()` using formula = 'y ~ x'

# save
ggsave(filename = "CEWL_temp_control_lm.pdf",
       plot = CEWL_control,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)
## `geom_smooth()` using formula = 'y ~ x'

CEWL (relative) ~ time

# n lizards
CEWL_time_series %>%
  dplyr::filter(complete.cases(time_min, relative_CEWL)) %>%
  group_by(lizard_ID) %>%
  summarise(max(relative_CEWL)) %>%
  nrow()
## [1] 37
# plot
ggplot() +
   geom_point(data = CEWL_time_series,
              aes(x= time_min, 
                  y = relative_CEWL,
                  shape = treatment,
                  color = treatment),
              size = 0.01, 
              alpha = 0.1) +
  
   geom_hline(yintercept = 0, size = 0.5,
              linetype = "dashed", color = "black") +
   
   geom_smooth(data = CEWL_time_series,
               aes(x = time_min, 
                   y = relative_CEWL,
                   color = treatment),
               method = "lm", 
               formula = y ~ poly((x), 2),
               size = 1.5, 
               se = F) +
   
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none") +
  
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  
  annotate(geom = "point", x = 15, y = -3.9, 
           size = 4, pch = 21, fill = control) +
   annotate(geom = "point", x = 15, y = -7.5, 
           size = 4, pch = 25, fill = cool) + 
   annotate(geom = "point", x = 15, y = 3, 
           size = 4, pch = 24, fill = heat) +
  scale_x_continuous(limits = c(1, 15), 
                      breaks = seq(1, 15, 2)) +

   labs(x = "Time (minutes)" , 
        y = bquote('CEWL (g '*m^-2*' '*h^-1*')')) -> CEWLr_time
  
CEWLr_time
## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 4 rows containing missing values (`geom_point()`).

# save
ggsave(filename = "CEWL_relative_time.pdf",
       plot = CEWLr_time,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 150, height = 100)
## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).
## Removed 4 rows containing missing values (`geom_point()`).

CEWL Rate of Change

nrow(lizard_coeffs)
## [1] 37
ggplot() +
  geom_jitter(data = lizard_coeffs,
              aes(x = treatment,
                   y = rate_change,
                   shape = treatment,
                   color = treatment,
                   fill = treatment),
                   size = 0.8,
               alpha = 0.3,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = emmeans,
                aes(x = treatment,
                    y = emmean, 
                    color = treatment,
                    group = treatment,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 0.9) +
  geom_point(data = emmeans, 
             aes(x = treatment,
                   y = emmean, 
                 shape = treatment,
                 fill = treatment), 
              color = "black",
              size = 4) +
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none") +
  
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
  scale_fill_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
  
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
  
   geom_hline(yintercept = 0, size = 0.5,
              linetype = "dashed", color = "black") +
  
   annotate("text", x = 1, y = 0.42, label = "A", size = 3) +
   annotate("text", x = 2, y = 0.3, label = "A", size = 3) +
   annotate("text", x = 3, y = 0.9, label = "B", size = 3) +
  
   labs(x = "" , 
        y = expression(Delta ~ 'CEWL '*min^-1*''),
        color = "Treatment") -> rate_change_CEWL_fig
rate_change_CEWL_fig

# save
ggsave(filename = "CEWL_change_min.pdf",
       plot = rate_change_CEWL_fig,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)

Experimental Design Skematic

exp_design <- read.csv("./Data/exp design.csv", 
                            header = TRUE,
                            fileEncoding="UTF-8-BOM")
ggplot(exp_design) +
   aes(x = time_min, 
       y = temp_C,
       color = treatment) +
     geom_line(size = 1) +
   geom_vline(xintercept = 60,
              alpha = 0.5,
              size = 0.5,
              linetype = "dashed", color = "black") +
   geom_text(aes(30, 37, family = "sans"), label = "a", color = "black", size = 3) +
   geom_vline(xintercept = 65,
              alpha = 0.5,
              size = 0.5,
              linetype = "dashed", color = "black") +
   geom_text(aes(62.5, 37, family = "sans"), label = "b", color = "black", size = 3) +
   geom_vline(xintercept = 80,
              alpha = 0.5,
              size = 0.5,
              linetype = "dashed", color = "black") +
   geom_text(aes(72.5, 37, family = "sans"), label = "c", color = "black", size = 3) +
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        axis.text.x = element_blank(),
        #axis.line.x = element_blank(),
        axis.ticks.x = element_blank()) +

   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time" , 
        y = "Body Temperature (°C)", 
        color = "Treatment") +
   
    annotate(geom = "point", x = 80, y = 25, 
           size = 3, pch = 21, fill = control) +
   annotate(geom = "point", x = 80, y = 20, 
           size = 3, pch = 25, fill = cool) + 
   annotate(geom = "point", x = 80, y = 30, 
           size = 3, pch = 24, fill = heat) +
   
   scale_y_continuous(limits = c(14, 37), 
                      breaks=seq(15, 35, 5)) -> methods_schmatic_fig
methods_schmatic_fig

# save
ggsave(filename = "methods_schematic.pdf",
       plot = methods_schmatic_fig,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)

Save Grouped Figures

Relative CEWL & Rate of Change

multi_fig_relative <- ggarrange(CEWLr_time, 
          ggarrange(rate_change_CEWL_fig, LEGEND, 
                    nrow = 1, ncol = 2,
                    widths = c(2,1),
                    align = "hv"
                    ),
          nrow = 2, ncol = 1,
          heights = c(2,1),
          labels = c("A", "B"))
## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 4 rows containing missing values (`geom_point()`).
multi_fig_relative

ggsave(filename = "CEWL_time_multi_fig.pdf",
       plot = multi_fig_relative,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 120, height = 160) # full width = 180

Raw CEWL ~ Cloacal temp & control zoom-in

multi_fig_CEWL_by_temp <- ggarrange(CEWL_ctemp, 
                  ggarrange(CEWL_control, LEGEND, 
                            nrow = 1, ncol = 2,
                            widths = c(2, 1),
                            align = "hv"
                            ),
                  nrow = 2, ncol = 1,
                  heights = c(2,1),
          labels = c("A", "B"))
## `geom_smooth()` using formula = 'y ~ x'
multi_fig_CEWL_by_temp

ggsave(filename = "CEWL_temp_multi_fig.pdf",
       plot = multi_fig_CEWL_by_temp,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 120, height = 160)